function y = runge_kutta4(a, b, n, y0,f)

	% Intrari:
		% a, b - intervalul de integrare
		% n - numarul de puncte
		% y0 - conditia initiala
		% f - functia y'=f(x,y)
	% Iesiri:
		% y - aproximarea solutiei in punctele xi

	h = ( b-a ) / n;
	y = zeros( n+1, 1 );
	y(1) = y0;
	for k = 2 : n+1
		x = a + k*h;
		K1 = feval(f, x, y(k-1));
		K2 = feval(f, x + h/2, y(k-1) + K1/2);
		K3 = feval(f, x + h/2, y(k-1) + K2/2);
		K4 = feval(f, x + h, y(k-1) + K3);
		y(k) = y(k-1)+h*(K1+2*K2+2*K3+K4)/6;
	endfor

endfunction
